# analysis.R

##### PACKAGES ####
# If you don't have these packages yet,
# install them using install.packages()
# before running the following library
# commands.

library(coefplot)
library(ggplot2)
library(tidyr)
library(ggpubr)
library(sjPlot)
library(interflex)
library(dplyr)
library(gt)
library(datawizard)
library(psych)
library(marginaleffects)

# Open the data ####
dat <- read.csv("CRT.csv")
dat_unchecked <- read.csv("CRT.csv")

# Preliminary variable recodes for descriptive statistics of sample ####
dat <- dat %>%
  mutate(gender_recode = 
           case_when(gender == 1 ~ "Male",
                     gender == 2 ~ "Female") %>%
           factor(levels = c("Male", "Female")),
         age_recode = 
           case_when(age >= 18 & age < 30 ~ "18-29",
                     age >= 30 & age < 45 ~ "30-44",
                     age >= 45 & age < 65 ~ "45-64",
                     age > 64 ~ "65 and over") %>%
           factor(levels = c("18-29", "30-44", "45-64", "65 and over")),
         race_recode = 
           case_when(ethnicity == 1 ~ "White",
                     ethnicity >= 2 & ethnicity <= 15 ~ "Non-white",
                     TRUE ~ "Would rather not answer") %>%
           factor(levels = c("White", "Non-white", "Would rather not answer")),
         education_recode = 
           case_when(education == 1 ~ "Some high school or less",
                     education == 2 ~ "High school graduate",
                     education == 3 ~ "Other post-high school vocational training",
                     education == 4 ~ "Completed some college, but no degree",
                     education == 5 ~ "Associate's degree",
                     education == 6 ~ "Bachelor's degree",
                     education == 7 ~ "Master's or professional degree",
                     education == 8 ~ "Doctorate degree",
                     education == -3105 ~ "None of the above  ") %>%
           factor(levels = c("Some high school or less",
                             "High school graduate",
                             "Other post-high school vocational training",
                             "Completed some college, but no degree",
                             "Associate's degree",
                             "Bachelor's degree",
                             "Master's or professional degree",
                             "Doctorate degree",
                             "None of the above  ")),
         party_recode = case_when(
           political_party == 1 | political_party == 2 |political_party == 3 |
             political_party == 6 ~ "Democrat",
           political_party == 10 | political_party == 9 | political_party == 8 |  
             political_party == 5 ~ "Republican",
           political_party == 4 ~ "Ind./Other",
           political_party == 7 ~ "Ind./Other") %>%
           factor(levels = c("Democrat", "Ind./Other","Republican")),
         ideology_recode = 
           case_when(political_party == 1 ~ "Strong Democrat",
                     political_party == 2 ~ "Not very strong Democrat",
                     political_party == 3 | political_party == 6 ~ "Leaning Democrat",
                     political_party == 4 ~ "Ind./Other ",
                     political_party == 5 | political_party == 8 ~ "Leaning Republican",
                     political_party == 9 ~ "Not very strong Republican",
                     political_party == 10 ~ "Strong Republican",
                     political_party == 7 ~ "Ind./Other") %>%
           factor(levels = c("Strong Democrat",
                             "Leaning Democrat",
                             "Not very strong Democrat",
                             "Ind./Other",
                             "Not very strong Republican",
                             "Leaning Republican",
                             "Strong Republican"))) 

# Create descriptive statistics table before drops (Table A5) ####
tableA5 <- 
  dat %>%
  tidyr::pivot_longer(cols = c(gender_recode,
                               age_recode,
                               race_recode,
                               education_recode,
                               party_recode,
                               ideology_recode)) %>%
  count(name, value, name = 'Percent') %>%
  group_by(name) %>%
  mutate(name = case_when(
    name == "age_recode" ~ "Age",
    name == "education_recode" ~ "Education",
    name == "gender_recode"~ "Gender",
    name == "ideology_recode" ~ "Ideology",
    name == "party_recode" ~ "Partisanship",
    name == "race_recode" ~ "Race")) %>%
  rename( Demographic = value) %>%
  mutate(Percent = prop.table(Percent)) %>% 
  gt() %>%
  tab_header("Table 3: Sample Pre-Attention Checks") %>%
  fmt_percent(columns = Percent) %>%
  tab_options(data_row.padding = px(1))


# Drop those who failed attention check and those who went too fast through the survey for main dataset ####
dat <- subset(dat, dat$fire._3==3)
dat <- subset(dat, dat$Duration..in.seconds.>150)

# Create descriptive statistics table after drops (Table A4) ####
tableA4 <- 
  dat %>%
  tidyr::pivot_longer(cols = c(gender_recode,
                               age_recode,
                               race_recode,
                               education_recode,
                               party_recode,
                               ideology_recode)) %>%
  count(name, value, name = 'Percent') %>%
  group_by(name) %>%
  mutate(name = case_when(
         name == "age_recode" ~ "Age",
         name == "education_recode" ~ "Education",
         name == "gender_recode"~ "Gender",
         name == "ideology_recode" ~ "Ideology",
         name == "party_recode" ~ "Partisanship",
         name == "race_recode" ~ "Race")) %>%
  rename( Demographic = value) %>%
  mutate(Percent = prop.table(Percent)) %>% 
  gt() %>%
  fmt_percent(columns = Percent) %>%
  tab_options(data_row.padding = px(1))


#### Topline results to pre-experiment CRT question -- statistics listed data/methods section text ####
prop.table(table(dat$crt_support))
prop.table(table(dat$crt_support, dat$party_recode),2)

##### Creating the Dependent Variables #####

# DV1 - vote
dat$localvotediff <- dat$vote_post_1 - dat$vote_pre._1
dat$statevotediff <- dat$vote_post_2 - dat$vote_pre._2
dat$combovotediff <- (dat$localvotediff + dat$statevotediff)/2
dat$meanvotepost <- (dat$vote_post_1 + dat$vote_post_2)/2
dat$meanvotepre <- (dat$vote_pre._1 + dat$vote_pre._1)/2

mean(dat$combovotediff, na.rm=T)

# DV2 - Encourage other to vote
dat$localdiff_encourage <- dat$encourage_post_1 - dat$encourage_pre_1
dat$statediff_encourage <- dat$encourage_post_2 - dat$encourage_pre_2
dat$combodiff_encourage <- (dat$localdiff_encourage + dat$statediff_encourage)/2
dat$meanencouragepost <- (dat$encourage_post_1 + dat$encourage_post_2)/2
dat$meanencouragepre <- (dat$encourage_pre_1 + dat$encourage_pre_2)/2

mean(dat$combodiff_encourage, na.rm=T)

# DV3 - Contact
dat$localdiff_contact <- dat$contact_post_1 - dat$contact_pre._1
dat$statediff_contact <- dat$contact_post_2 - dat$contact_pre._2
dat$combodiff_contact <- (dat$localdiff_contact + dat$statediff_contact)/2
dat$meancontactpost <- (dat$contact_post_1 + dat$contact_post_2)/2
dat$meancontactpre <- (dat$contact_pre._1 + dat$contact_pre._2)/2

mean(dat$combodiff_contact, na.rm=T)

# DV4 - Campaign/volunteer
dat$localdiff_campaign <- dat$volunteer_post_1 - dat$campaign_pre._1
dat$statediff_campaign <- dat$volunteer_post_2 - dat$campaign_pre._2
dat$combodiff_campaign <- (dat$localdiff_campaign + dat$statediff_campaign)/2
dat$meancampaignpost <- (dat$volunteer_post_1 + dat$volunteer_post_2)/2
dat$meancampaignpre <- (dat$campaign_pre._1 + dat$campaign_pre._2)/2

mean(dat$combodiff_campaign, na.rm=T)

# DV5 - Run
dat$localdiff_run <- dat$run_post_1 - dat$run_pre_1
dat$statediff_run <- dat$run_post_2 - dat$run_pre_2
dat$combodiff_run <- (dat$localdiff_run + dat$statediff_run)/2
dat$meanrunpost <- (dat$run_post_1 + dat$run_post_2)/2
dat$meanrunpre <- (dat$run_pre_1 + dat$run_pre_2)/2

mean(dat$combodiff_run, na.rm=T)

# DV6 - Index (create index for all pre and post variables using factor analysis)
fitpre <- fa(dat[c("vote_pre._1", "vote_pre._2", "encourage_pre_1", "encourage_pre_2", "contact_pre._1", "contact_pre._2", "campaign_pre._1", "campaign_pre._2", "run_pre_1", "run_pre_2")], fm="pa")
fitpre
dat$indexpre <- predict(fitpre, dat[c("vote_pre._1", "vote_pre._2", "encourage_pre_1", "encourage_pre_2", "contact_pre._1", "contact_pre._2", "campaign_pre._1", "campaign_pre._2", "run_pre_1", "run_pre_2")])
dat$indexpre <- datawizard::rescale(dat$indexpre, to = c(0, 100))

fitpost <- fa(dat[c("vote_post_1", "vote_post_2", "encourage_post_1", "encourage_post_2", "contact_post_1", "contact_post_2", "volunteer_post_1", "volunteer_post_2", "run_post_1", "run_post_2")], fm="pa")
fitpost
dat$indexpost <- predict(fitpost, dat[c("vote_post_1", "vote_post_2", "encourage_post_1", "encourage_post_2", "contact_post_1", "contact_post_2", "volunteer_post_1", "volunteer_post_2", "run_post_1", "run_post_2")])
dat$indexpost <- datawizard::rescale(dat$indexpost, to = c(0, 100))

# Create table of factor loadings for Supplementary Information - Table A1 ####
loadings <- as.data.frame(cbind(fitpre$loadings, fitpost$loadings))
loadings <- setNames(loadings,c("Pre-treatment","Post-treatment"))
loadings$items <- c("Vote-local", "Vote-state", "Encourage-local", "Encourage-state", "Contact-local", "Contact-state", "Campaign-local", "Campaign-state", "Run-local", "Run-state")
loadings

# Create new treatment indicators ####
dat$venue2[dat$venue=="city or town council"] <- "Local"
dat$venue2[dat$venue=="state legislature"] <- "State"
dat$venue2 <- factor(dat$venue2, levels=c("State", "Local"))

dat$action <- NA
dat$action[dat$outcome=="decided to mandate that public schools include critical race theory in their curricula"] <- "Teach CRT"
dat$action[dat$outcome=="decided to prohibit public schools from including critical race theory in their curricula"] <- "Ban CRT"
dat$action[dat$outcome=="decided that they would not take any action on the issue of critical race theory"] <- "No action"
dat$action <- factor(dat$action, levels=c("No action", "Teach CRT", "Ban CRT"))


# Creating factor variable for support ####
# Excluding not sure
dat$crt_support2 <- NA
dat$crt_support2[dat$crt_support==1] <- "Support"
dat$crt_support2[dat$crt_support==2] <- "Oppose"
dat$crt_support2 <- factor(dat$crt_support2, c("Support", "Oppose"))

# Including not sure
dat$crt_support3 <- NA
dat$crt_support3[dat$crt_support==1] <- "Support"
dat$crt_support3[dat$crt_support==2] <- "Oppose"
dat$crt_support3[dat$crt_support==3] <- "Not sure"
dat$crt_support3 <- factor(dat$crt_support3, c("Support", "Oppose", "Not sure"))

# Figure 4 - support/opposition for action taken in vignette ####

all <- as.data.frame(round(prop.table(table(dat$action, dat$treatment),1),2)*100)
all$group <- "All respondents"
reps <- as.data.frame(round(prop.table(table(dat$action[dat$party_recode=="Republican"], dat$treatment[dat$party_recode=="Republican"]),1),2)*100)
reps$group <- "Republicans"
dems <- as.data.frame(round(prop.table(table(dat$action[dat$party_recode=="Democrat"], dat$treatment[dat$party_recode=="Democrat"]),1),2)*100)
dems$group <- "Democrats"

party <- rbind(reps, dems)
party$group <- factor(party$group, levels=c("Republicans", "Democrats"))

allplot <- ggplot(all, aes(x=Var2, y=Freq, fill=group)) + 
  geom_bar(position="dodge", stat="identity") + ylim(0,60) +
  scale_x_discrete(labels=c("Strongly\nsupport", "Somewhat\nsupport", "Somewhat\noppose", "Strongly\nOppose", "Not\nsure")) +
  scale_fill_brewer(palette="Dark2") + facet_wrap(~Var1) +
  ggtitle("All respondents") + xlab("") + ylab("Percent") +
  theme_minimal() + theme(legend.position="none")

partyplot <- ggplot(party, aes(x=Var2, y=Freq, fill=group)) + 
  geom_bar(position="dodge", stat="identity") + ylim(0,60) +
  scale_x_discrete(labels=c("Strongly\nsupport", "Somewhat\nsupport", "Somewhat\noppose", "Strongly\nOppose", "Not\nsure")) +
  scale_fill_brewer(palette="Set1") + facet_wrap(~Var1) +
  ggtitle("By partisanship") + xlab("") + ylab("Percent") +
  theme_minimal() + theme(legend.position="bottom", legend.title=element_blank())


png(filename="fig4.png", res=100, width=960, height=600)
ggarrange(allplot, partyplot, ncol=1)
dev.off()

# Figure 5: Treatment effect of hypothetical government action on CRT on likelihood of taking action, conditional on opinion towards CRT ####

# Vote
vote <- lm(meanvotepost ~ action*crt_support3 + meanvotepre, data=dat)
votep <- marginaleffects(vote, by="crt_support3", variables="action")

voteplot <- ggplot(votep, aes(y=estimate, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nworking for campaign") + 
  ggtitle("Vote") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-11, 22)

# Encouragement
encourage <- lm(meanencouragepost ~ action*crt_support3 + meanencouragepre, data=dat)
encouragep <- marginaleffects(encourage, by="crt_support3", variables="action")

encourageplot <- ggplot(encouragep, aes(y=estimate, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nworking for campaign") + 
  ggtitle("Encourage") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-11, 22)

# Contact
contact <- lm(meancontactpost ~ action*crt_support3 + meancontactpre, data=dat)
contactp <- marginaleffects(contact, by="crt_support3", variables="action")

contactplot <- ggplot(contactp, aes(y=estimate, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nworking for campaign") + 
  ggtitle("Contact") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-11, 22)

# Campaign
campaign <- lm(meancampaignpost ~ action*crt_support3 + meancampaignpre, data=dat)
campaignp <- marginaleffects(campaign, by="crt_support3", variables="action")

campaignplot <- ggplot(campaignp, aes(y=estimate, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nworking for campaign") + 
  ggtitle("Campaign") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-11, 22)

# Run
run <- lm(meanrunpost ~ action*crt_support3 + meanrunpre, data=dat)
runp <- marginaleffects(run, by="crt_support3", variables="action")

runplot <- ggplot(runp, aes(y=estimate, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nrunning for office") + 
  ggtitle("Run") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-11, 22)

# Index
index <- lm(indexpost ~ action*crt_support3 + indexpre, data=dat)
indexp <- marginaleffects(index, by="crt_support3", variables="action")

indexplot <- ggplot(indexp, aes(y=estimate, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on\nparticipation index") + 
  ggtitle("Participation index") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-11, 22)

# Pull individual grobs together and make final graph
combogrid <- ggarrange(voteplot, encourageplot, contactplot, campaignplot, runplot, indexplot, common.legend = TRUE) + theme_minimal()

combogrid

png(filename="fig5.png", res=100, width=960, height=600)
combogrid
dev.off()

# Holm p-value correction
alltests <- rbind(votep %>% mutate(model = "vote"), 
                  encouragep %>% mutate(model = "encourage"), 
                  contactp %>% mutate(model = "contact"), 
                  campaignp %>% mutate(model = "campaign"), 
                  runp %>% mutate(model = "run") # , 
                  #indexp # index is a collection of DVs not independent of each other
) 

alltests$adjusted.p.values <- round(p.adjust(alltests$p.value, method="holm"),4) 

alltests <- alltests %>%
  select(model, contrast, crt_support3, p.value, adjusted.p.values) %>%
  arrange(adjusted.p.values)


# Table 1: OLS regressions testing effect of experimental conditions on likelihood of engaging in political activities
 ####

m1 <- lm(vote_post_1 ~ action + vote_pre._1, data=dat)
m2 <- lm(vote_post_2 ~ action + vote_pre._2, data=dat)
m3 <- lm(volunteer_post_1 ~ action + campaign_pre._1, data=dat)
m4 <- lm(volunteer_post_2 ~ action + campaign_pre._2, data=dat)
m5 <- lm(run_post_1 ~ action + run_pre_1, data=dat)
m6 <- lm(run_post_2 ~ action + run_pre_2, data=dat)
m7 <- lm(contact_post_1 ~ action + contact_pre._1, data=dat)
m8 <- lm(contact_post_2 ~ action + contact_pre._2, data=dat)
m9 <- lm(encourage_post_1 ~ action + encourage_pre_1, data=dat)
m10 <- lm(encourage_post_2 ~ action + encourage_pre_2, data=dat)

tab_model(
  m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, 
  dv.labels = c("Vote\nlocal", "Vote\nstate", "Campaign\nlocal", "Campaign\nstate", "Run\nlocal", "Run\nstate", "Contact\nlocal", "Contact\nstate", "Encourage\nlocal", "Encourage\nstate"),
  pred.labels = c("Intercept", "Teach CRT", "Ban CRT"),
  show.ci = FALSE, 
  show.p = FALSE, 
  show.se = TRUE,
  collapse.se = TRUE,
  p.style = c("stars"),
  p.threshold = c(0.05),
  transform = NULL,
  file = "table1.html"
)

# Holm p-value correction
alltests_table1 <- rbind(summary(m1)$coefficients[,4], 
                         summary(m2)$coefficients[,4], 
                         summary(m3)$coefficients[,4], 
                         summary(m4)$coefficients[,4], 
                         summary(m5)$coefficients[,4], 
                         summary(m6)$coefficients[,4], 
                         summary(m7)$coefficients[,4], 
                         summary(m8)$coefficients[,4], 
                         summary(m9)$coefficients[,4], 
                         summary(m10)$coefficients[,4]) %>%
  as.data.frame() %>%
  mutate(model = c("vote-local",
                   "vote-state",
                   "campaign-local",
                   "campaign-state",
                   "run-local",
                   "run-state",
                   "contact-local",
                   "contact-state",
                   "encourage-local",
                   "encourage-state")) %>% 
  pivot_longer(cols = !(model)) %>%
  filter(name != "(Intercept)" & name!= "vote_pre._1") # Intercepts are not p-values we pay attention
                                # to in our hypotheses
  
alltests_table1$adjusted.p.values <- round(p.adjust(alltests_table1$value, method="holm"),4)

alltests_table1 <- alltests_table1 %>% rename(p.value = value)
# write.csv(alltests_table1, "alltests_table1.csv")

# Table A2: OLS regression testing effect of experimental conditions on likelihood of engaging in political activities including interactions for venue ####

m1 <- lm(vote_post_1 ~ action*venue2 + vote_pre._1, data=dat)
m2 <- lm(vote_post_2 ~ action*venue2 + vote_pre._2, data=dat)
m3 <- lm(volunteer_post_1 ~ action*venue2 + campaign_pre._1, data=dat)
m4 <- lm(volunteer_post_2 ~ action*venue2 + campaign_pre._2, data=dat)
m5 <- lm(run_post_1 ~ action*venue2 + run_pre_1, data=dat)
m6 <- lm(run_post_2 ~ action*venue2 + run_pre_2, data=dat)
m7 <- lm(contact_post_1 ~ action*venue2 + contact_pre._1, data=dat)
m8 <- lm(contact_post_2 ~ action*venue2 + contact_pre._2, data=dat)
m9 <- lm(encourage_post_1 ~ action*venue2 + encourage_pre_1, data=dat)
m10 <- lm(encourage_post_2 ~ action*venue2 + encourage_pre_2, data=dat)

tab_model(
  m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, 
  dv.labels = c("Vote\nlocal", "Vote\nstate", "Campaign\nlocal", "Campaign\nstate", "Run\nlocal", "Run\nstate", "Contact\nlocal", "Contact\nstate", "Encourage\nlocal", "Encourage\nstate"),
  pred.labels = c("Intercept", "Teach CRT", "Ban CRT", "Venue = local", "Teach x local", "Ban x local"),
  show.ci = FALSE, 
  show.p = FALSE, 
  show.se = TRUE,
  collapse.se = TRUE,
  p.style = c("stars"),
  p.threshold = c(0.05),
  transform = NULL,
  file = "tableA2.html"
)


# Figure A3: Reproduction of Figure 2 separating by local and state action ####

# Local
# Vote
vote <- lm(vote_post_1 ~ action*crt_support3 + vote_pre._1, data=dat)
votep <- marginaleffects(vote, by="crt_support3", variables="action")

voteplot <- ggplot(votep, aes(y=estimate, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nvoting") + 
  ggtitle("Vote") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-15, 22)

# Encouragement
encourage <- lm(encourage_post_1 ~ action*crt_support3 + encourage_pre_1, data=dat)
encouragep <- marginaleffects(encourage, by="crt_support3", variables="action")

encourageplot <- ggplot(encouragep, aes(y=estimate, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nencouraging others to vote") + 
  ggtitle("Encourage") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-15, 22)

# Contact
contact <- lm(contact_post_1 ~ action*crt_support3 + contact_pre._1, data=dat)
contactp <- marginaleffects(contact, by="crt_support3", variables="action")

contactplot <- ggplot(contactp, aes(y=estimate, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\ncontacting elected officials") + 
  ggtitle("Contact") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-15, 22)

# Campaign
campaign <- lm(volunteer_post_1 ~ action*crt_support3 + campaign_pre._1, data=dat)
campaignp <- marginaleffects(campaign, by="crt_support3", variables="action")

campaignplot <- ggplot(campaignp, aes(y=estimate, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nworking for campaign") + 
  ggtitle("Campaign") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-15, 22)

# Run
run <- lm(run_post_1 ~ action*crt_support3 + run_pre_1, data=dat)
runp <- marginaleffects(run, by="crt_support3", variables="action")

runplot <- ggplot(runp, aes(y=estimate, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nrunning for office") + 
  ggtitle("Run") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-15, 22)


# Pull individual grobs together
combogrid <- ggarrange(voteplot, encourageplot, contactplot, campaignplot, runplot, common.legend = TRUE) + theme_minimal() + ggtitle("Local Participation Measures")

png(filename="figA3panel1.png", res=100, width=960, height=600)
combogrid
dev.off()

# State
# Vote
vote <- lm(vote_post_2 ~ action*crt_support3 + vote_pre._2, data=dat)
votep <- marginaleffects(vote, by="crt_support3", variables="action")

voteplot <- ggplot(votep, aes(y=dydx, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nvoting") + 
  ggtitle("Vote") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-15, 25)

# Encouragement
encourage <- lm(encourage_post_2 ~ action*crt_support3 + encourage_pre_2, data=dat)
encouragep <- marginaleffects(encourage, by="crt_support3", variables="action")

encourageplot <- ggplot(encouragep, aes(y=dydx, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nencouraging others to vote") + 
  ggtitle("Encourage") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-15, 25)

# Contact
contact <- lm(contact_post_2 ~ action*crt_support3 + contact_pre._2, data=dat)
contactp <- marginaleffects(contact, by="crt_support3", variables="action")

contactplot <- ggplot(contactp, aes(y=dydx, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\ncontacting elected officials") + 
  ggtitle("Contact") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-15, 25)

# Campaign
campaign <- lm(volunteer_post_2 ~ action*crt_support3 + campaign_pre._2, data=dat)
campaignp <- marginaleffects(campaign, by="crt_support3", variables="action")

campaignplot <- ggplot(campaignp, aes(y=dydx, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nworking for campaign") + 
  ggtitle("Campaign") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-15, 25)

# Run
run <- lm(run_post_2 ~ action*crt_support3 + run_pre_2, data=dat)
runp <- marginaleffects(run, by="crt_support3", variables="action")

runplot <- ggplot(runp, aes(y=dydx, x=crt_support3, group=contrast, shape=contrast)) + 
  geom_point(stat="identity", position=position_dodge(width=.2)) + 
  geom_errorbar(aes(x=crt_support3, ymin=conf.low, ymax=conf.high), position=position_dodge(width=.2), width=.2) +
  xlab("Subject's position on CRT") +
  ylab("Treatment effect on likelihood of\nrunning for office") + 
  ggtitle("Run") + theme_minimal() + 
  scale_shape_manual('Treatment effect of:',values=c(1, 2),labels=c("Banning CRT", "Teaching CRT")) +
  theme(legend.position="bottom") +
  geom_hline(yintercept=0) + ylim(-15, 25)


# Pull individual grobs together
combogrid <- ggarrange(voteplot, encourageplot, contactplot, campaignplot, runplot, common.legend = TRUE) + theme_minimal() + ggtitle("State Participation Measures")

combogrid

png(filename="figA3panel3.png", res=100, width=960, height=600)
combogrid
dev.off()

# Figures 6 and 7: treatment effects by partisanship and racism denial ####

# Create racism denial scale 
dat$denial1 <- car::recode(dat$fire._1, "1=0; 2=.2; 3=.4; 4=.6; 5=.8; 6=1")
dat$denial2 <- car::recode(dat$fire._2, "1=1; 2=.8; 3=.6; 4=.4; 5=.2; 6=0")
dat$denial <- rowMeans(dat[c("denial1", "denial2")])

# Create treatment variable that is zero for control and 1 for teach CRT 
dat$teachCRT <- NA
dat$teachCRT[dat$action=="No action"] <- 0
dat$teachCRT[dat$action=="Teach CRT"] <- 1

# Run models
contact <- lm(meancontactpost ~  teachCRT*denial + teachCRT*party_recode + meancontactpre  , dat)
vote <- lm(meanvotepost ~  teachCRT*denial + teachCRT*party_recode + meanvotepre, dat)
campaign <- lm(meancampaignpost ~  teachCRT*denial + teachCRT*party_recode + meancampaignpre, dat)
encourage <- lm(meanencouragepost ~  teachCRT*denial + teachCRT*party_recode + meanencouragepre , dat)
run <- lm(meanrunpost ~  teachCRT*denial + teachCRT*party_recode + meanrunpre, dat)
index <- lm(indexpost ~ teachCRT*denial + teachCRT*party_recode + indexpre, dat)


# Holm p-value correction
alltests_racismdenial<- rbind(summary(contact)$coefficients[,4], 
                        summary(campaign)$coefficients[,4], 
                        summary(encourage)$coefficients[,4], 
                        summary(vote)$coefficients[,4],
                        summary(run)$coefficients[,4] #,
                        # summary(index)$coefficients[,4]
                        ) %>%
  as.data.frame() %>%
  mutate(model = c("contact", "campaign", "encourage", "vote", 
                   "run" #, 
                   #"index"
                   )) %>%
  pivot_longer(cols = !(model)) %>%
  filter(name != "(Intercept)" & name != "meancontactpre" & 
           name != "meanvotepre" & name != "meanencouragepre" & 
           name != "meanrunpre" & name != "meancampaignpre")
        # removing lagged controls

alltests_racismdenial$adjusted.p.values <- round(p.adjust(alltests_racismdenial$value, method="holm"),4)

# write.csv(alltests_racismdenial, "alltests_racismdenial.csv")

# Create Figure 7 ####
contactp <- plot_cme(contact, variables = "teachCRT", condition = list("denial", "party_recode"="Republican")) + 
  ggtitle("Contact") + ylab("Treatment effect of\n teaching CRT") + 
  scale_color_manual(values=c("black")) + scale_fill_manual(values=c("black")) +
  xlab("Denial of racism") + geom_hline(yintercept = 0) + ylim(-20, 35) + 
  theme_minimal() + theme(legend.position="none") +
  scale_x_continuous(breaks=c(0, .5,  1)) 

campaignp <- plot_cme(campaign, variables = "teachCRT", condition = list("denial", "party_recode"="Republican")) + 
  ggtitle("Campaign") + ylab("Treatment effect of\n teaching CRT") +
  scale_color_manual(values=c("black")) + scale_fill_manual(values=c("black")) +
  xlab("Denial of racism") + geom_hline(yintercept = 0) + ylim(-20, 35) +
  scale_x_continuous(breaks=c(0, .5,  1)) + 
  theme_minimal() + theme(legend.position="none") 

encouragep <- plot_cme(encourage, variables = "teachCRT", condition = list("denial", "party_recode"="Republican")) + 
  ggtitle("Encourage") + ylab("Treatment effect of\n teaching CRT") +
  scale_color_manual(values=c("black")) + scale_fill_manual(values=c("black")) +
  xlab("Denial of racism") + geom_hline(yintercept = 0) + ylim(-20, 35) +
  scale_x_continuous(breaks=c(0, .5,  1)) + 
  theme_minimal() + theme(legend.position="none") 

votep <- plot_cme(vote, variables = "teachCRT", condition = list("denial", "party_recode"="Republican")) + 
  ggtitle("Vote") + ylab("Treatment effect of\n teaching CRT") +
  scale_color_manual(values=c("black")) + scale_fill_manual(values=c("black")) +
  xlab("Denial of racism") + geom_hline(yintercept = 0) + ylim(-20, 35)  +
  scale_x_continuous(breaks=c(0, .5,  1)) + 
  theme_minimal() + theme(legend.position="none") 

runp <- plot_cme(run, variables = "teachCRT", condition = list("denial", "party_recode"="Republican")) + 
  ggtitle("Run") + ylab("Treatment effect of\n teaching CRT") +
  scale_color_manual(values=c("black")) + scale_fill_manual(values=c("black")) +
  xlab("Denial of racism") + geom_hline(yintercept = 0) + ylim(-20, 35)  +
  scale_x_continuous(breaks=c(0, .5,  1)) + 
  theme_minimal() + theme(legend.position="none") 

indexp <- plot_cme(index, variables = "teachCRT", condition = list("denial", "party_recode"="Republican")) + 
  ggtitle("Index") + ylab("Treatment effect of\n teaching CRT") +
  scale_color_manual(values=c("black")) + scale_fill_manual(values=c("black")) +
  xlab("Denial of racism") + geom_hline(yintercept = 0) + ylim(-20, 35)  +
  scale_x_continuous(breaks=c(0, .5,  1)) + 
  theme_minimal() + theme(legend.position="none") 

racismgrid <- ggarrange(votep, encouragep, contactp, campaignp, runp, indexp) + theme_minimal()

png(filename="fig7.png", res=100, width=960, height=600)
racismgrid
dev.off()

# Create Figure 6 ####

contactp2 <- plot_cme(contact, variables = "teachCRT", condition = c("party_recode")) + 
  ggtitle("Contact") + ylab("Treatment effect of\n teaching CRT") +
  geom_hline(yintercept = 0) + ylim(-10, 20) + xlab("") + 
  scale_x_discrete(labels=c("Dem.", "Ind.", "Rep.")) + theme_minimal()

campaignp2 <- plot_cme(campaign, variables = "teachCRT", condition = c("party_recode")) + 
  ggtitle("Campaign") + ylab("Treatment effect of\n teaching CRT") +
  geom_hline(yintercept = 0) + ylim(-10, 20) + xlab("") +
  scale_x_discrete(labels=c("Dem.", "Ind.", "Rep.")) + theme_minimal()

encouragep2 <- plot_cme(encourage, variables = "teachCRT", condition = c("party_recode")) + 
  ggtitle("Encourage") + ylab("Treatment effect of\n teaching CRT") +
  geom_hline(yintercept = 0) + ylim(-10, 20) + xlab("") +
  scale_x_discrete(labels=c("Dem.", "Ind.", "Rep.")) + theme_minimal()

votep2 <- plot_cme(vote, variables = "teachCRT", condition = c("party_recode")) + 
  ggtitle("Vote") + ylab("Treatment effect of\n teaching CRT") +
  geom_hline(yintercept = 0) + ylim(-10, 20) + xlab("") +
  scale_x_discrete(labels=c("Dem.", "Ind.", "Rep.")) + theme_minimal()

runp2 <- plot_cme(run, variables = "teachCRT", condition = c("party_recode")) + 
  ggtitle("Run") + ylab("Treatment effect of\n teaching CRT") +
  geom_hline(yintercept = 0) + ylim(-10, 20) + xlab("") +
  scale_x_discrete(labels=c("Dem.", "Ind.", "Rep.")) + theme_minimal()

indexp2 <- plot_cme(index, variables = "teachCRT", condition = c("party_recode")) + 
  ggtitle("Index") + ylab("Treatment effect of\n teaching CRT") +
   geom_hline(yintercept = 0) + ylim(-10, 20) + xlab("") +
  scale_x_discrete(labels=c("Dem.", "Ind.", "Rep.")) + theme_minimal()


partygrid <- ggarrange(votep2, encouragep2, contactp2, campaignp2, runp2, indexp2) 

png(filename="fig6.png", res=100, width=960, height=600)
partygrid
dev.off()


# Figure A4: Interflex plots for SI ####
contact <- interflex(Y = "meancontactpost", D = "teachCRT", X = "denial", Z="meancontactpre", data = dat, estimator = "binning", na.rm=T, nbins=3)
contact <- plot(contact, ylab="Treatment effect", xlab = "Denial of racism", theme.bw = T, main = "Contact", ylim=c(-20, 40)) 
vote <- interflex(Y = "meanvotepost", D = "teachCRT", X = "denial", Z="meanvotepre", data = dat, estimator = "binning", na.rm=T, nbins=3)
vote <- plot(vote, ylab="Treatment effect", xlab = "Denial of racism", theme.bw = T, main = "Vote", ylim=c(-20, 40)) 
campaign <- interflex(Y = "meancampaignpost", D = "teachCRT", X = "denial", Z="meancampaignpre", data = dat, estimator = "binning", na.rm=T, nbins=3)
campaign <- plot(campaign, ylab="Treatment effect", xlab = "Denial of racism", theme.bw = T, main = "Campaign", ylim=c(-20, 40)) 
encourage <- interflex(Y = "meanencouragepost", D = "teachCRT", X = "denial", Z="meanencouragepre", data = dat, estimator = "binning", na.rm=T, nbins=3)
encourage <- plot(encourage, ylab="Treatment effect", xlab = "Denial of racism", theme.bw = T, main = "Encourage", ylim=c(-20, 40)) 
run <- interflex(Y = "meanrunpost", D = "teachCRT", X = "denial", Z="meanrunpre", data = dat, estimator = "binning", na.rm=T, nbins=3)
run <- plot(run, ylab="Treatment effect", xlab = "Denial of racism", theme.bw = T, main = "Run", ylim=c(-20, 40)) 
index <- interflex(Y = "indexpost", D = "teachCRT", X = "denial", Z="indexpre", data = dat, estimator = "binning", na.rm=T, nbins=3)
index <- plot(index, ylab="Treatment effect", xlab = "Denial of racism", theme.bw = T, main = "Index", ylim=c(-20, 40)) 

png(filename="figA4.png", res=100, width=960, height=600)
ggarrange(vote, encourage, contact, campaign, run, index, common.legend = TRUE) + theme_minimal() + theme(legend.title=element_blank())
dev.off()


